home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Precision Software Appli…tions Silver Collection 4
/
Precision Software Applications Silver Collection Volume 4 (1993).iso
/
stats
/
chadyn.exe
/
YDBLRTR.C
< prev
next >
Wrap
C/C++ Source or Header
|
1988-11-28
|
5KB
|
163 lines
/********************* YDblRtr.C for the double rotor ************************/
/********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
#include "yinclud.h"
double exp();
static double u11,
u12,
u21,
u22,
s1,
s2;
/* for double rotor: */
static double Veloc11,
Veloc12,
Veloc21,
Veloc22;
static double Matrix11,
Matrix12,
Matrix21,
Matrix22;
static double cc1,
cc2;
start33 () { /* calculate parameters for the double rotor
map */
double T,
mass1,
mass2,
length1,
length2,
damping1,
damping2;
/* above parameters mostly to be set from menu */
double inertia1,
inertia2;
double discrim,
expon1,
expon2,
unorm;
mass1 = C1; /* mass1 */
mass2 = C2; /* mass2 */
length2 = C4; /* length2 -- length1 is set below */
damping1 = C5; /* damping1 & damping2 are damping */
damping2 = C6;
T = C7;
/* Let the length of the first rod and the inertias be set to values that
allow for a symmetric set of ODEs for the system */
length1 = length2 * sqrt(mass2 / (mass1 + mass2));
inertia1 = (mass1 + mass2) * length1 * length1;
inertia2 = mass2 * length2 * length2;
/* compute eigenvalues for revised map */
discrim = 0.5 * sqrt(damping1 * damping1 + 4.0 * damping2 * damping2);
s1 = 0.5 * damping1 + damping2 + discrim;
s2 = 0.5 * damping1 + damping2 - discrim;
u11 = damping2 - s1;
u12 = damping2;
unorm = sqrt(u11 * u11 + u12 * u12);
u11 = u11 / unorm;
u12 = u12 / unorm;
u21 = damping2 - s2;
u22 = damping2;
unorm = sqrt(u21 * u21 + u22 * u22);
u21 = u21 / unorm;
u22 = u22 / unorm;
/* Compute matrix M(T) for theta dot */
expon1 = exp(-s1 * T);
expon2 = exp(-s2 * T);
Veloc11 = u11 * u11 * expon1 + u21 * u21 * expon2;
Veloc12 = u11 * u12 * expon1 + u21 * u22 * expon2;
Veloc21 = u11 * u12 * expon1 + u21 * u22 * expon2;
Veloc22 = u12 * u12 * expon1 + u22 * u22 * expon2;
/* Compute matrix M(T) for theta */
Matrix11 = -(u11 * u11 * (expon1 - 1) / s1 + u21 * u21 * (expon2 - 1) / s2);
Matrix12 = -(u11 * u12 * (expon1 - 1) / s1 + u21 * u22 * (expon2 - 1) / s2);
Matrix21 = -(u11 * u12 * (expon1 - 1) / s1 + u21 * u22 * (expon2 - 1) / s2);
Matrix22 = -(u12 * u12 * (expon1 - 1) / s1 + u22 * u22 * (expon2 - 1) / s2);
/* miscellaneous constants -- rho is forcing and is only in iterated map */
cc1 = length1 / inertia1;
cc2 = length2 / inertia2;
}
dblrotor() {
double moduloAB();
int i;
int base;
double y_temp[20],
yval0, /* y0, y1 are global */
yval1;
yval0 = y[0] + Matrix11 * y[2] + Matrix12 * y[3];
y_temp[0] = yval0;
yval1 = y[1] + Matrix21 * y[2] + Matrix22 * y[3];
y_temp[1] = yval1;
y_temp[2] = Veloc11 * y[2] + Veloc12 * y[3] + rho * cc1 * sin(yval0);
y_temp[3] = Veloc21 * y[2] + Veloc22 * y[3] + rho * cc2 * sin(yval1);
for(i = 0; i < num_lyap; i++) {
base = lyapzero + vec_dim * i;
y_temp[base]
= y[base] + Matrix11 * y[base + 2] + Matrix12 * y[base + 3];
y_temp[base + 1]
= y[base + 1] + Matrix21 * y[base + 2] + Matrix22 * y[base + 3];
y_temp[base + 2]
= Veloc11 * y[base + 2] + Veloc12 * y[base + 3]
+ rho * cc1 * cos(yval0) * (y[base]
+ Matrix11 * y[base + 2] + Matrix12 * y[base + 3]);
y_temp[base + 3]
= Veloc21 * y[base + 2] + Veloc22 * y[base + 3]
+ rho * cc2 * cos(yval1) * (y[base + 1] +
Matrix21 * y[base + 2] + Matrix22 * y[base + 3]);
}
for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
y[i] = y_temp[i];
/* now bring y[0] and y[1] into -pi,pi etc. */
y[0] = moduloAB(y[0], -pi, pi);
y[1] = moduloAB(y[1], -pi, pi);
}
initDblRotor() {
extern int start33 ();
extern int dblrotor();
vec_dim = 4; /* the dimension of the Lyapunov vectors =
phase space dim */
num_lyap = 0; /* the number of Lyapunov numbers to be
computed <= vec_dim */
lyapzero = 4; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
init_process = start33;
dim = lyapzero + num_lyap * vec_dim;
/* needed for rungekutta */
X_upper = pi; /* x scale */
X_lower = -pi;
Y_upper = pi; /* y scale */
Y_lower = -pi;
/* mass1,2 = C1,C2; length2 = C4; length1=C4*sqrt(C2/ (C1+C2));
damping1,2 = C5,C6 hit time=C7 */
C1 = 1.;
C2 = 1.;
C4 = 1.;
C5 = 1.;
C6 = 1.;
C7 = 1.;
rho = 9.; /* forcing strength */
map = dblrotor;
}